Preface

The goal of this assignment is to develop some intuition about the impact of the number of nodes in the hidden layer of the neural network. We will use few simulated examples to have clear understanding of the structure of the data we are modeling and will assess how performance of the neural network model is impacted by the structure in the data and the setup of the network.

First of all, to compensate for lack of coverage on this topic in ISLR, let’s go over a couple of simple examples. We start with simulating a simple two class dataset in 2D predictor space with an outcome representative of an interaction between attributes.

# fix seed so that narrative always matches the plots:
set.seed(1234567890)
nObs <- 1000
ctrPos <- 2
xyTmp <- matrix(rnorm(4*nObs),ncol=2)
xyCtrsTmp <- matrix(sample(c(-1,1)*ctrPos,nObs*4,replace=TRUE),ncol=2)
xyTmp <- xyTmp + xyCtrsTmp
gTmp <- paste0("class",(1+sign(apply(xyCtrsTmp,1,prod)))/2)
plot(xyTmp,col=as.numeric(factor(gTmp)),pch=as.numeric(factor(gTmp)),xlab="X1",ylab="X2")
abline(h=0)
abline(v=0)

Symbol color and shape indicate class. Typical problem that will present a problem for ant approach estimating a single linear decision boundary. We used similar simulated data for the random forest assignment.

One hidden node

We can fit simple neural network (using all default values in the call to neuralnet – notice that both covariates and outcome have to be numeric as opposed to factor) and plot its layout (allowing for its output to be included in Rmarkdown generated report actually seems to be quite painful - one has to overwrite original implementation of plot.nn with the one that doesn’t call dev.new() that is included in this Rmarkdown file with echo=FALSE – to do the same you have to include that block into your Rmarkdown file also):

### Doesn't run: "requires numeric/complex ... arguments"
### nnRes <- neuralnet(g~X1+X2,data.frame(g=gTmp,xyTmp))
nnRes <- neuralnet(g~X1+X2,data.frame(g=as.numeric(factor(gTmp)),xyTmp))
plot(nnRes)

That shows us model one node in a single hidden layer (default parameters).

We can lookup actual model predictions and recalculate them from input variables (in the field covariate) and model weight and activation function (fields weights and act.fct respectively):

head(nnRes$net.result[[1]])
##          [,1]
## 1 1.632380797
## 2 1.632380951
## 3 1.632376796
## 4 1.547455068
## 5 1.632380951
## 6 1.632380931
cbind(rep(1,6),nnRes$act.fct(cbind(rep(1,6),nnRes$covariate[1:6,])%*%nnRes$weights[[1]][[1]]))%*%nnRes$weights[[1]][[2]]
##             [,1]
## [1,] 1.632380797
## [2,] 1.632380951
## [3,] 1.632376796
## [4,] 1.547455068
## [5,] 1.632380951
## [6,] 1.632380931

Notice that input parameter linear.output governs whether activation function is called on the value of the output node or not:

nnResNLO <- neuralnet(g~X1+X2,data.frame(g=as.numeric(factor(gTmp)),xyTmp),linear.output=FALSE)
head(nnResNLO$net.result[[1]])
##           [,1]
## 1 0.9999903688
## 2 0.9999903531
## 3 0.9999903914
## 4 0.9999904043
## 5 0.9999903216
## 6 0.9999903639
nnResNLO$act.fct(cbind(rep(1,6),nnResNLO$act.fct(cbind(rep(1,6),nnResNLO$covariate[1:6,])%*%nnResNLO$weights[[1]][[1]]))%*%nnResNLO$weights[[1]][[2]])
##              [,1]
## [1,] 0.9999903688
## [2,] 0.9999903531
## [3,] 0.9999903914
## [4,] 0.9999904043
## [5,] 0.9999903216
## [6,] 0.9999903639

We can plot model output indicating class identity (left panel below) that tells us that majority of the points were estimated to be close to 1 or 1.6 and that majority of those estimated to be close to 1 correspond to the about half of the observations at the first level of the class factor (i.e. numerical value of 1). It is also easy to see that those with predicted value of 1.6 represent roughly 1:2 mix of observations from the first and second levels of the outcome respectively, so that \(1.6 \approx (1+2*2)/3 = 5/3\) approximately equals average of their numerical values corresponding to the levels of the factor representing them.

The nature of the model estimated by neuralnet in this (very simple!) case becomes even more intuitive if we render all points in the area encompassing our training set with model predictions and overlay training dataset on top of that (right panel below). It is immediately apparent that this model identified a line in this 2D space separating one cloud of points belonging mostly to one class from all others so that predicted values are approximately equal to the average outcome on each side of this decision boundary.

plotNNpreds2D2class <- function(inpNN,inpClassThresh,inpGrid=(-60:60)/10) {
  plot(inpNN$net.result[[1]],col=inpNN$response,pch=inpNN$response)
  table(inpNN$net.result[[1]][,1]>inpClassThresh,inpNN$response)
  xyGridTmp <- cbind(X1=rep(inpGrid,length(inpGrid)),X2=sort(rep(inpGrid,length(inpGrid))))
  gridValsTmp <- compute(inpNN,xyGridTmp)
  plot(xyGridTmp,col=as.numeric(gridValsTmp$net.result>inpClassThresh)+1,pch=20,cex=0.3,main=paste("Error:",signif(sum(inpNN$err.fct(inpNN$response,inpNN$net.result[[1]][,1])),6)))
  points(inpNN$covariate,col=inpNN$response,pch=inpNN$response)
  ## Equations defining decision boundary:
  ## 1*w0 + X1*w1 + X2*w2 = 0, i.e.:
  ## 0 = inpNN$weights[[1]][1]+inpNN$weights[[1]][2]*X1+inpNN$weights[[1]][3]*X2, i.e:
  ## X2 = (-inpNN$weights[[1]][1] - inpNN$weights[[1]][2]*X1) / inpNN$weights[[1]][3]
  for ( iTmp in 1:ncol(inpNN$weights[[1]][[1]]) ) {
    abline(-inpNN$weights[[1]][[1]][1,iTmp] / inpNN$weights[[1]][[1]][3,iTmp], -inpNN$weights[[1]][[1]][2,iTmp] /inpNN$weights[[1]][[1]][3,iTmp],lty=2,lwd=2)
  }
}
old.par <- par(mfrow=c(1,2),ps=16)
plotNNpreds2D2class(nnRes,1.3)

par(old.par)

The important points resulting from the figure above are the following:

  • as simple of a model as the one that was employed here (with one node in a single hidden layer along with all other default parameters) used cannot do much better than what we observed here
  • because calling (default - logistic) activation function on a given linear combination of the input variables more or less amounts to assigning almost all points on one side of hyperplane (line in 2D, plane in 3D, etc.) to zero and on the other side – to unity
  • weights involved in transforming outcome of the hidden layer into model predictions will change those zeroes and ones to values closer to the desired outcome values, but still, use if such a simple model (with a single hidden node) employed here to prime our intuition more or amounts to spliting covariate space into two half spaces by a hyperplane and assigning constant outcomes to the vast majority of the points on either side of it
  • the weights for the inputs into the single hidden node shown in the network layout plot above and stored in weights field of the result returned by neuralnet define this hyperplane (line in 2D, etc.) shown as dashes in the right panel above
  • this hyperplane is where sum of weighted input variables and a constant is identical zero (and thus the result of logistic activation function is 0.5 rapidly becoming zero or one for points further away from this boundary)

Two hidden nodes, single hidden layer

Now, let’s add another node to the hidden layer of this network. From the above, we know what to expect as a result of that – another hyperplane (line in 2D) will be added to the space of covariates now dividing it into (depending on whether those hyperplans are almost parallel or not) three or four subspaces, consequently, assigning most of the points to three or four potentially different constants. Clearly, this level of granularity could suffice for developing a model that would do quite well in our toy example.

To have more than one node in the hidden layer we set hidden parameter to the number of nodes in it (length of vector provided as hidden parameter governs the number of the hidden layers in the network – we still use one layer here):

nnRes2 <- neuralnet(g~X1+X2,data.frame(g=as.numeric(factor(gTmp)),xyTmp),hidden=2)
plot(nnRes2)

We can see that now resulting network has two nodes in a single hidden layer, which two covariates enter with weights that are approximately comparable in magnitude and opposite in sign. Their comparably weighted sum added to a constant close to one gives the outcome value of this model. The effect of those weights in defining decision boundaries in the space of predictors is best seen from the figure below:

old.par <- par(mfrow=c(1,2),ps=16)
plotNNpreds2D2class(nnRes2,1.5)

par(old.par)

This model sets up two almost parallel lines that encompass most of the observations from the first class, leaving most of the observations from the second class outside of the resulting slab. Now let’s repeat fitting neural network three times (each time starting with random choice of starting weights in the model) and compare stability of the resulting models:

old.par <- par(mfcol=c(2,3),ps=16)
for ( iTry in 1:3 ) {
  nnRes2 <- neuralnet(g~X1+X2,data.frame(g=as.numeric(factor(gTmp)),xyTmp),hidden=2)
  plotNNpreds2D2class(nnRes2,1.5)
}

par(old.par)

We can see that quite frequently given the parameters used the process converges to a suboptimal solution with about half of the observations remaining in “gray” zone where their assignment to either of the classes is not immediately apparent.

Aside from the multitude of local minima for neural network fitting procedure that could prevent it from finding better solutions, the main point to take from this exercise is that adding more nodes to the hidden layer (with all other default choices employed here) amounts to adding more hyperplanes bisecting the space of predictors, creating more and more subspaces where the outcome can take different values (often close to a constant in each subspace). Obviously, the geometry of the resulting decision surfaces can become quite complicated even with modest number of nodes in the hidden layer. Lastly, these considerations provide some intuition for considering what could be a useful number of hidden nodes in the model. In thinking about that it might be useful to consider how many such hyperplanes could be sufficient to effectively separate observations belonging to different outcome categories. Not that we necessarily would have such knowledge ahead of time, but this might prove to be a complementary way to think about the problem in addition to the often sited empirical guidelines that are based on the number of predictor variables, etc.

The point of this homework is to assess how these aspects of neural network fitting play out in another simulated dataset.

Problem 1 (10 points): 3D data with spherical class boundary

Simulate data with n=1000 observations and p=3 covariates – all random variables from standard normal distribution. Create two category class variable assigning all observations within a sphere with radius of 1.5 centered at 3D zero to one class category and all others – to the second. Since you will be reusing this code in the following two problems it is probably best to turn this procedure into a function with appropriate parameters. Check that resulting class assignments split these observations very roughly evenly between these two groups. Plot values of the resulting covariates projected at each pair of the axes indicating classes to which observations belong with symbol color and/or shape (you can use function pairs, for example). What is the smallest number of planes in 3D space that would completely enclose points from the “inner” class?

simData <- function(n=1000, cols=3, inform=3, radius=1.5){
  set.seed(1234)
  x <- matrix(rnorm(n*cols),ncol=cols)
  cl<-as.numeric(factor(sqrt(rowSums(x[,1:inform]^2))<radius))
  #points within inner circle are set to class 2, outside =1  
  data.frame(CLASS=cl, x)
}


d1 <- simData()
head(d1)
##   CLASS            X1            X2             X3
## 1     1 -1.2070657494 -1.2053334196 -0.97381861936
## 2     2  0.2774292421  0.3014667398 -0.09963119709
## 3     1  1.0844411767 -1.5391452236 -0.11073500011
## 4     1 -2.3456977026  0.6353707235  1.19219459601
## 5     1  0.4291246888  0.7029517748 -1.65588593019
## 6     1  0.5060558922 -1.9058828646 -1.04564329068
pairs(d1[2:4], main = "", pch=21,bg = c("red","green3")[d1$CLASS])

#3d Visualisation
attach(d1)
sp <- scatterplot3d(X1,X2,X3,color=CLASS,cex.symbols=0.8, pch=19,angle=20)

#4 planes would be the minimum required to completely enclose a sphere

Problem 2 (20 points): neural network classifier

For the dataset simulated above fit neural networks with 1 through 6 nodes in a single hidden layer (use neuralnet implementation). For each of them calculate training error (see an example in Preface where it was calculated using err.fct field in the result returned by neuralnet). Simulate another independent dataset (with n=10,000 observations to make resulting test error estimates less variable) using the same procedure as above (3D, two classes, decision boundary as a sphere of 1.5 radius) and use it to calculate test error at each number of hidden nodes. Plot training and test errors as function of the number of nodes in the hidden layer. What does resulting plot tells you about the interplay between model error, model complexity and problem geometry? What is the geometrical interpretation of this error behavior?

#Create test dataset
test <- simData(n=10000)

train.errs <- numeric(6)
test.errs  <- numeric(6)
for (i in 1:6){
  nn <- neuralnet(CLASS~X1+X2+X3,data=d1,err.fct="sse",hidden=i)
  #plot(nn)
  #Given the difference in observations between test and train, we will conver to MSE
  train.errs[i] <- nn$result.matrix[1]/nrow(d1)
  my.predict <- compute(nn, test[,2:4])$net.result
  test.errs[i] <- sum((test$CLASS-my.predict)^2)/nrow(test)
}

#Reshape data for ggplot
dd <- data.frame(rbind(train.errs,test.errs))
dd$lbl <- rownames(dd)
dx <-gather(dd,idx,err,X1:X6)
ggplot(dx,aes(idx,err,group=lbl)) + geom_line(aes(linetype=lbl, color=lbl)) + geom_point(aes(color=lbl)) + labs(x="Hidden Nodes" , y="MSE") +
  theme(legend.position="top")

#We can see that the MSE drops until 4 nodes, then we see a levelling off of MSE despite an increase in model complexity. Geometrically, this is what we assumed in problem 1: That we would need a minimum of 4 hyperplanes to enclose our classification boundary that is a sphere.

Problem 3 (30 points): evaluate impacts of sample size and noise

Setup a simulation repeating procedure described above for n=100, 200 and 500 observations in the training set as well adding none, 1, 2 and 5 null variables to the training and test data (and to the covariates in formula provided to neuralnet). Draw values for null variables from standard normal distribution as well and do not use them in the assignment of the observations to the class category (e.g. x<-matrix(rnorm(600),ncol=6); cl<-as.numeric(factor(sqrt(rowSums(x[,1:3]^2))<1.5)) creates dataset with three informative and three null variables). Repeat calculation of training and test errors at least several times for each combination of sample size, number of null variables and size of the hidden layer simulating new training and test dataset every time to assess variability in those estimates. Present resulting error rates so that the effects of sample size and fraction of null variables can be discerned and discuss their impact of the resulting model fits.

dfTmp <- NULL
for (ns in c(100,200,500)){
  for (nulls in c(0,1,2,5)){
    pop <- simData(n=ns, cols=3+nulls)
    nm <- names(pop)
    fm <- as.formula(paste("CLASS ~", paste(nm[!nm %in% "CLASS"], collapse = " + ")))
    for (hidden in 1:6){
    
        bTrain <- sample(c(FALSE,TRUE),nrow(pop),replace=TRUE)
        train <- pop[bTrain,]
        test <- pop[!bTrain,]

        nn <- neuralnet(fm,data=train,err.fct="sse",hidden=hidden)
        predict <- compute(nn, test[,-1])$net.result

        dfTmp <- rbind(dfTmp,data.frame(train.err=nn$result.matrix[1]/nrow(train),
                   test.err=sum((test$CLASS-predict)^2)/nrow(test),
                   hnodes=hidden,
                   nulls=nulls,
                   n=ns))
      
    }
  }
}

dx <- gather(dfTmp, err.type, err.val, train.err:test.err, factor_key=TRUE)

ggplot(dx,aes(x=as.factor(hnodes),y=err.val,shape=err.type,colour=n))+geom_point()+facet_wrap(~nulls) 

#We can see that an increase in the number of covariates has a large effect on the test error when n is small. As sample size increases, test error starts decreasing. Training error is harldy affected by this.

#

Extra 20 points problem: model banknote authentication data

Use neuralnet to model the outcome in banknote authentication dataset that we used in previous assignments and compare its test error at several sizes of hidden layer to that observed for SVM and KNN approaches.

dbaDat <- read.table("data_banknote_authentication.txt",sep=",")
colnames(dbaDat) <- c("var","skew","curt","entr","auth")

cls <- as.factor(dbaDat$auth)
set.seed(1234)

#Using 4 nodes per hidden layer, as 4 nodes performed the best previously
#10 recplicates to see variance of errors
dfTmp <- NULL
for (i in 1:10){
  for (j in 1:10){
    bTrain <- sample(c(FALSE,TRUE),nrow(dbaDat),replace=TRUE)
    train <- dbaDat[bTrain,]
    test <- dbaDat[!bTrain,]
    nn <- neuralnet(auth~var+skew+curt+entr,err.fct="sse",data=train,hidden=rep(4,i))
    cm <- as.matrix(table(test$auth,1+(compute(nn, test[,1:4])$net.result>1.5)))
    dfTmp <- rbind(dfTmp, data.frame(i,err= 1-(sum(diag(cm))/sum(cm))))
  }
}
boxplot(err~i, data=dfTmp, xlab="Hidden Layers")

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:neuralnet':
## 
##     compute
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
dfTmp %>% group_by(i) %>% summarize(mean=mean(err),sd=sd(err))
## # A tibble: 10 × 3
##        i         mean             sd
##    <int>        <dbl>          <dbl>
## 1      1 0.4405856692 0.015311929671
## 2      2 0.4471207983 0.019908265364
## 3      3 0.4425642358 0.009544485724
## 4      4 0.4465815719 0.016111468991
## 5      5 0.4497809521 0.011835899853
## 6      6 0.4417826854 0.011443185908
## 7      7 0.4441066346 0.017088008294
## 8      8 0.4513915247 0.014306502543
## 9      9 0.4542380697 0.010237581837
## 10    10 0.4431740796 0.011621295049
#dplyr compute has interactions with neuralnet, so detaching
detach("package:dplyr", unload=TRUE)
#We can see that the best performance for 4 nodes is given with 2 hidden layers

set.seed(1234)
bTrain <- sample(c(FALSE,TRUE),nrow(dbaDat),replace=TRUE)
train <- dbaDat[bTrain,]
test <- dbaDat[!bTrain,]

#Test misclassifcation rates
testError <- function(truth, predicted){
  cm <- as.matrix(table(truth,predicted))
  1-(sum(diag(cm))/sum(cm))
}



nm <- neuralnet(as.numeric(factor(auth))~var+skew+curt+entr,err.fct="sse",data=train,hidden=rep(4,2))
nn.err <- testError(test$auth,
                    1+(compute(nm, test[,1:4])$net.result>1.5))
plot(nm)

km <- knn(train,test,k=2,cl=cls[bTrain])
kn.err <- testError(test$auth,km)

svm <- svm(as.factor(auth)~var+skew+curt+entr, data=train, kernel="polynomial", gamma=.2,
cost =20)

svm.err <- testError(test$auth,predict(svm
                                       ,test[,-5]))


#Testing against best performing KNN and SVM models from Week 11
library(knitr)
kable(data.frame(nn.err, kn.err, svm.err))
nn.err kn.err svm.err
0 0 0.0102639296 
#The performance of the neuralnet seems highly variable to the initial seed provided, similar to KNN. However, on this test set nn and knn performed equally well with SVM performing slightly worse